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Abstract 

We investigate the effects of disorder on single particle time-evolution and two-particle corre- 
lations in an array of evanescently coupled waveguides with position-dependent tunneling rates. 
In the clean limit, the energy spectrum of such an array is widely tunable. In the presence of a 
Hermitian on-site or tunneling disorder, we find that the localization of a wave packet is highly 
sensitive to this energy spectrum. In particular, for an input confined to a single waveguide, we 
show that the fraction of light localized to the original waveguide depends on the tunneling profile. 
We compare the two-particle intensity correlations in the presence of Hermitian, tunneling disorder 
and non-Hermitian, parity-and-time-reversal {VT) symmetric, on-site potential disorder. We show 
the two-particle correlation function in both cases is qualitatively similar, since both disorders 
preserve the particle-hole symmetric nature of the energy spectrum. 



I. INTRODUCTION 



The development of integrated photonic structures in recent years has allowed researchers 
to study many condensed matter effects in optical systems. For example, the motion of 
an electron in a periodic potential has been translated into the propagation of light in a 
waveguide array and Bloch-oscillation-like phenomena have been directly observed ^IH3]. 
Optical waveguide arrays have also attracted considerable attention because the diffraction 
properties of light in these structures are very different from that in bulk media |1] . 

In recent years, arrays of evanescently coupled waveguides have attracted considerable 
interest as versatile structures suitable for the study of quantum and condensed matter phe- 
nomena such as Bloch oscillations, Zener tunneling [5], Dirac zitterbewegung [6J, quantum 
random walks [7J, and quantum ratchets [8]. This interest is a result of being able to ob- 
serve quantum phenomena on a length-scale of a few millimeters, the typical length of a 
waveguide, as the discrete Schrodinger equation for a particle on a lattice is identical to the 
longitudinal component of the Helmholtz equation for the electric field in an array of waveg- 
uides P, [in] • Non-classical states of light have also been used to study quantum properties 
of light propagating through waveguide arrays [2 [3l [H [TT] . In particular, the quantum cor- 
relations of two non-interacting indistinguishable particles propagating simultaneously have 
shown the effect of their initial separation on propagation of the particles in clean |T2] and 
disordered [13] systems. 

Over the last decade, since the seminal work of Bender and co-workers p3l [15] . there has 
been much interest in the study of systems that are described by non-Hermitian, but parity- 
{V) and time-reversal- (T) symmetric Hamiltonians. This symmetry has led to predictions 
of new phenomenon such as Bloch oscillations in complex crystals [16], an optical medium 
that can simultaneously act as an emitter and a perfect absorber of coherent waves [17] . 
and induced quantum coherence between Bose- Einstein condensates [18]. PT-symmetry 
breaking in a classical system has recently been experimentally observed in waveguide arrays 
[T^EU]. Some of the predictions for waveguide arrays with PT-symmetric Hamiltonians are 
Bloch oscillations in arrays with defects [2T] , the invisibility of defects in such an array [22] , 
directed transport in non-linear arrays [23], periodic wave packet reconstruction [21], and 
double refraction [25] . Most of the waveguide arrays considered in these works are uniform, 
meaning the evanescent couphng between adjacent waveguides is the same. However, over 
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the past year, a strong interest in waveguide arrays with position-dependent tunnehng rate 
or evanescent couphng has developed [211 [23, EZ]; this interest is driven, in part, by the 
tremendous flexibihty that is offered by optical waveguides where the evanescent coupling 
can be varied by two orders of magnitude, and by the tunable properties of wave packet 
evolution in such arrays with non-uniform, position-dependent tunneling rates |27j . 

Recent theoretical and experimental reports have explored the role of disorder in cou- 
pled waveguide arrays, and demonstrated Anderson localization of light in one-dimension. 
However, there are still several fundamental questions that remain unresolved with regard 
to the localization phenomenon (see [28]). For example, how is the localization defined in a 
lattice which departs from the standard configuration in which the tunneling amplitudes are 
equal? By studying different models of non-uniform tunneling, one can probe the nature of 
localization. 

Finally, note that although there is an exact mapping between coupled optical waveguides 
and a one-dimensional lattice model, qualitatively, coupled optical waveguides allow explo- 
ration of dynamics and disorder effects in a parameter regime that is virtually inaccessible 
in their condensed matter counterparts: the number of waveguides ~ 10 — 100 is much 
smaller than typical number of lattice sites ^ 10^ in any quantum wire; the typical disor- 
der strength A in optical waveguides is comparable to the bandwidth, whereas in condensed 
matter systems, A < Ep ^ bandwidth; the typical wavepacket in coupled waveguides spans 
the entire bandwidth whereas in condensed matter systems, a typical wavepacket only spans 
energies comparable to Ep <^ bandwidth. Thus, it is crucially important to quantitatively 
define and explore the interaction between dynamics and disorder in lattice models with a 
position-dependent tunneling amplitude that are applicable to optical waveguide lattices. 

The inclusion of position-dependent tunneling rates, and the presence of disorder, makes 
our problem analytically intractable. Therefore, most studies on disorder effects on light 
propagation in waveguide arrays are numerical in nature. The few exceptions are the works 
on Bloch oscillations and related effects, in which the disorder is absent, and the tunneling 
amplitudes are constant. Since the fabrication of waveguides with non- uniform tunneling 
rates is technologically feasible, as is the inclusion of disorder in such arrays, it is important 
to investigate the evolution of light in such systems. As we show in this paper, the results 
are dramatically different from what one finds in systems with constant tunneling. 

In this paper, we investigate the effect of disorder in waveguide arrays with a position- 
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dependent tunneling rate C{j). We study the time evolution of a wave packet and the effect 
the tunneling rate's functional form has upon the disorder-induced localization of the wave 
packet. The appropriate functional form of the tunneling function increases the localization 
due to disorder, but, for a large disorder strength, the fraction of light localized to the initial 
waveguide may be independent of the global structure of the array. We find that, based on 
the initial wave function profile, disorder can cause broadening or localize the wave packet to 
two different waveguides. We compare and contrast the intensity-intensity correlations for 
random, Hermitian, off-diagonal disorder and random, non-Hermitian, PT-symmetric, loss 
and gain disorder. We show that although both disorders give rise to qualitatively identical 
correlation functions, the full correlation matrix can distinguish between the two. 

The plan of the paper is as follows. In Sec. |Tl] we introduce the tight-binding lattice 



Hamiltonian, the tunneling functions, and the disorder we consider. Sec. Ill contains results 
for the competition between tunneling and disorder, and how it affects the time-evolution and 



localization of a wave packet. In Sec. IV we extend the model to include non-Hermitian, 
PT-symmetric, loss and gain disorder, and present the results for the intensity-intensity 
correlations that arise from two distinct disorders. We conclude the paper with a brief 
discussion in Sec. |Vl 



II. TIGHT-BINDING MODEL 

We consider an array of single-mode, evanescently coupled waveguides where the prop- 
agation of light is described by a tight-binding Hamiltonian. In second quantized form, the 
Hamiltonian is given by 

N N-1 

where Oj {aj) represents the creation (annihilation) operator for a photon in waveguide j, 
f3j is the linear-propagation constant (or equivalently the potential) at site j, and C{j) is 
the tunneling rate between waveguides j and j + 1. A uniform array is characterized by 
a constant tunneling rate, C{j) = Cq. We choose a parity-symmetric tunneling function 
characterized by a single parameter a [27] 

Co.ij) = Co[jiN-j)r/\ (2) 
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Experimentally, one can engineer this tunneling function by symmetrically increasing (de- 
creasing) the center-to-center distance between adjacent waveguides for negative (positive) 
values of a. When a > 0, the tunneling rate is maximum at the center of the array 
whereas for a < 0, it is minimum. The energy-spectrum bandwidth is defined as 
Aq, = -Emax — -Emin whcrc Ej are eigenvalues of the Hamiltonian ([T|; for a clean system, it 
scales as A«(A^) ~ ^CqA^" for a > and A^, ~ /iCqA^-I^I/^ for a < [26]. We choose the 
bandwidth inverse as the characteristic time-scale for the system, Ta{N) = h/Aa{N). 

We introduce a Hermitian disorder through random variations of the on-site potential 
The mean value of on-site potential is irrelevant as long as it is the same for each waveguide; 
it only introduces a constant shift of the energy spectrum, and does not affect the bandwidth 
or the characteristic time. Therefore, it is set to zero. Since the disorder-induced localization 
is independent of the probability distribution of disorder provided the different distributions 
have zero mean and the same variance we use a Gaussian distribution, 

where is the variance of the distribution, and thus, ha characterizes the strength of the 
disorder. Note that, in contrast with the waveguide-dependent tunneling rate, the variance 
of disorder is not dependent upon the waveguide index j. 

Starting with an arbitrary initial state |V'(0)) and the Hermitian Hamiltonian Ha, Eq.([l]), 
with given on-site disorder potentials, we obtain the unitary time-evolution operator U (t) = 
exp{—iHt/h), the time-evolved wave function lipit)) = U{t)\ilj{0)) , and the time-evolved, 
waveguide- index dependent intensity I{j,t) = |(j|-?/'(t))p. We consider an array with ^ 1 
waveguides, and use Box-MuUer algorithm to generate the random disorder with zero mean 
and the desired variance [22] • The characteristic time and the energy bandwidth A^ 
change with each disorder realization, and disorder averaging is carried out over A^^ ^ 
10^ — 10^ realizations to ensure that the results are independent of the number of realizations. 



III. WAVEPACKET EVOLUTION IN THE PRESENCE OF DISORDER 

In a uniform waveguide array, a weak disorder exponentially localizes a wave packet that 
is initially confined to a single waveguide [9]. We begin by investigating the effect of position- 
dependent tunneling function Ca{j) on the localization of light that is initially confined to 
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a single waveguide by calculating the a-dependence of the fraction of total intensity that 
remains in the initial waveguide at times t ^ for varying strengths of disorder. Similar to 
the energy-scale, disorder ha is measured in units of clean-system bandwidth Aa \ although, 
we will see below that the question of whether to use a global-scale, such as the bandwidth, 
or a local-scale, such as the tunneling rate, is essentially determined by the position and the 
profile of the initial wave packet. 
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FIG. 1. (Color Online) The left panel shows the steady-state intensity localized to the initial 
waveguide as a function of a for an array with = 100 waveguides, A',. = 10® disorder realizations, 
and input state \ip{0)) = \mo). The input state is given by mo = 50 (blue circles) and mo = 15 (red 
squares), and the disorder strength is higher than the bandwidth ha/ Aa^ = 3. We see that the 
localized fraction is weakly dependent upon a. The right panel shows corresponding results for an 
array with = 37 waveguides, A',. = 10^ disorder realizations, and a weaker disorder ha/ Aa^ = 1. 
We see that the localized fraction as a function of a depends acutely on different initial input 
states, mo = 11 (blue squares) and mo = 5 (red squares), when the input locations are relatively 
close to the boundary. 

The left panel in Fig. [l] shows the fraction of total intensity that remains in the initial 
waveguide as a function of a for = 100, initial waveguides mo = 50 (blue circles) and 
mo = 15 (red squares), and disorder strength ha/Aa^ = 3. When the initial waveguide is 
near the center of the array, for a < 0, we see that the localized fraction rapidly saturates 
as |a| increases, whereas for a > 0, the localized fraction is approximately independent of 
a. On the other hand, when the initial waveguide is near the edge, mo = 15, we see a clear 
dependence of the localized fraction on the tunneling exponent a. Note that for a weak 
disorder, ha/Aa^ <^ 1, the saturation of the localized fraction occurs for a < 0, whereas for 
a > the disorder has minimal effect and the average intensity at a particular waveguide is 
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nearly uniform, given by /a = The right panel in Fig. [T] shows the localized fraction as 
a function of a for a smaller disorder strength ha/Aa^ = 1, with = 37 and different input 
waveguide locations mo = 11 (blue circles) and ttiq = 5 (red squares). In general, we see that 
the a-dependence of the steady-state localized fraction is sensitive to the proximity of the 
initial waveguide location to the edge of the array. This result shows that for non-uniform 
waveguide arrays, the disorder-averaged steady-state intensity profile for a given input state 
|mo) is primarily determined by the competition between disorder and the local tunneling 
rate C(mo), instead of disorder and the clean-system bandwidth Ai°\ When ha ^ Aa \ we 
find that the localized intensity fraction in the initial waveguide is independent of a since the 
ballistic component of the intensity time-evolution is completely suppressed by the disorder. 
Finally, we note that the results displayed in Fig. [T] have been systematically checked by us 
to confirm that the effects shown are not due to the size of the lattice, and are due to the 
relative proximity of the input wavepacket to the boundary. 

We now consider a-dependence of the disorder-averaged intensity profile I{j, t) for a weak 
disorder fia/Aa^ = 0.05 when the input state is localized in waveguide mo = 15 near the 
edge of an = 100 waveguide array, \ip{0)) = |mo). Figure |2] shows that initially, the wave 
packet ballistically spreads and then, at times t/Ta > 300 ^ 1, develops a steady-state 
intensity profile that reflects its localization. The top panel (a) shows that for a uniform 
array, a = 0, the intensity is maximum in the initial waveguide, and decays monotonically 
away from it; note that the maximum intensity /(mo) ~ 0.02 is approximately twice the 
average intensity = 1/A^ = 0.01. The middle panel (b) shows that for a = 1, the 
disorder- averaged intensity profile changes from periodic reconstruction to localization as 
time increases. In contrast with the uniform tunneling a = case, the steady-state intensity 
in this case shows two peaks, with different weights, at mirror-symmetric positions 
mo = 15 and + 1 — mo = 86. The bottom panel (c) shows that for a = 2, the intensity 
profile has two sharply defined peaks at mirror symmetric locations. We emphasize that 
the existence of two peaks in the disorder-averaged intensity is a generic feature of parity- 
symmetric tunneling function Ca, and the ratio of weights at the two peaks can be varied 
by changing the disorder strength. Thus, in contrast with the exponential localization in a 
uniform waveguide array, the localization profile of a wave packet in an array with a ^ can 
be varied dramatically by appropriate choice of the input waveguide mo and disorder strength 
ha. When the disorder strength is increased further, ha/ Aa^ > 1, the twin-peak structure 
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FIG. 2. (Color Online) a-dependent disorder averaged intensity I{j,t) for an input state \ip{0)) = 
\mQ = 15) in an array with N = 100 waveguides, a weak disorder ha/ Aa^ = 0.05, and Nr = 10^ 
disorder realizations. The horizontal axis in each panel indicates time in units of Tq. Panel (a) 
shows exponential localization with a single peak at mo for an array with uniform tunneling a = 0. 
Panels (b) and (c) show corresponding results for a = 1 and a = 2 respectively. In each case, 
a ballistic expansion and reconstruction is followed by emergence of steady-state intensity profile 
that has two peaks, one at the input waveguide mo and the other at its mirror-symmetric 
counterpart, + 1 — mo- The relative weights at the two peaks can be tuned by the varying the 
weak disorder and the input position mo- 

disappears and the localization profile approaches an exponential as is expected [TT] . 

The results so far have been on the evolution of a single particle that is input to a single 
waveguide. In recent years, motivated primarily by the desire to understand the quantum 
correlations between indistinguishable particles, there have been studies on the evolution of 
two particles in waveguide lattices. Bromberg and co-workers have investigated quantum 
and classical correlations when two photons are coupled to either the same waveguide or to 
adjacent waveguides |3], and Lahini and co-workers [12] have studied a similar problem in 
disordered lattices. More recently, coupling of two phase-displaced wavepackets into adjacent 
waveguides has been used to propose a form of the quantum ratchet [8] . Motivated by these 
works, we next explore the effects of disorder and a-dependent clean-system spectrum on 
an input wave packet that is localized to two waveguides with a relative phase 6 between 
them, \ijjg{0)) = (\p) + e'^|g)) /\/2, with I < p,q < N. The phase 9 determines the time- 
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evolved intensity I{j,t) in a clean system and for a = 1, the phase information can only 
be extracted from intensity measurements in certain time-windows [2S] • Figure [3] shows the 
interplay between the phase 6 and weak disorder, and their effect on the disorder-averaged 
steady state intensity profile for a non- uniform waveguide array with = 60, a = 1, and 
an initial state with p = 20 and q = N — p = 40. 




normalized time waveguide number 

FIG. 3. (Color Online) Left-hand panel shows intensity I{j,t) as a function of 9 for an array with 
= 60 waveguides, disorder ha/ Aa^ = 0.05, and Nj- = 10^ disorder realizations. The initial input 
state is \ip{0)) = (|20) + e*^|40))/-v/2. The top, middle, and bottom panels correspond to = 0, 
6 = Tr/2 and 6 = tt respectively. The interference structure at short times t/ra ^ 100 is replaced 
by a steady-state double-peak intensity profile at times t/Ta> 300. The peak intensity is twice the 
average intensity la = ^/N = 0.0167 per waveguide. Right-hand panel shows the corresponding 
steady-state intensity at t/ra = 600 as a function of phase 6. The intensity near the waveguide 
array center shows enhancement for 9 = (blue solid line) and suppression for 9 = ir (black dotted 
line) when compared with the corresponding intensity for = 7r/2 (red dashed line). As expected, 
this intensity difference vanishes for a strong disorder. 

The left-hand panel shows the disorder-averaged intensity I{j,t) for = (top panel), 
6 = 7r/2 (middle panel), and 9 = tt (bottom panel) with a weak disorder, ha/Aa '^ = 0.05, 
and A'',. = 10^. At short times t/ra ^ 300, the intensity shows clear signatures of 6'-dependent 
interference and reconstruction due to equidistant energy levels of a clean a = 1 system. At 
large times t/Ta > 300, a steady-state intensity profile with two equal-weight, broad peaks 
near waveguides p = 20 and q = N — p = 40 emerges. A visual inspection of the three 
intensity profiles shows that the phase information is encoded in the steady-state intensity 
near the center of the waveguide array. The right-hand panel shows the steady-state intensity 
profile at time t/ra = 600 as a function of the phase 9] recaU that the average intensity 
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is given by la = ^/N = 0.0167. We see clearly that ^ = (blue solid line) and 6 = n (black 
dotted line) are marked by increased and suppressed intensity at the center of the waveguide 
array respectively, compared to the intensity value for 6 = 7i/2 (red dashed line). Thus, the 
phase information, accessible only in certain time windows in clean system, can be extracted 
from the disorder-averaged steady-state intensity profile; as the disorder gets stronger, this 
phase information disappears. 

These results show that a complete characterization of the disorder-induced steady-state 
intensity profile and the time at which steady-state is achieved is a nontrivial problem. 
In particular, we find that the steady-state intensity I{j) depends acutely on the size and 
relative phase 6 of the initial input state {ipelO)), the tunneling function a, and the proximity 
of the input state with the boundaries. 

IV. INTENSITY CORRELATIONS IN THE PRESENCE OF DISORDER: HER- 
MITIAN VS. NON-HERMITIAN CASE 

In the previous section, we focused on localization due to a Gaussian, Hermitian disorder 
in the on-site potential because the localization intensity profile is independent of the disorder 
origin (on-site potential or tunneling rates) and disorder probability distribution as long as 
different distributions have zero mean and the same variance [HI |30]. In contrast to this, 
higher-order intensity correlations provide a clue into the origin of the disorder; this is 
because an on-site disorder destroys the particle-hole symmetry of the clean-system energy 
spectrum whereas a tunneling-rate disorder preserves it [3jDJ. In this section, we explore 
the effects of a non-Hermitian, PT-symmetric, on-site disorder on the intensity-intensity 
correlations. To this end, we modify the Hamiltonian Ha to include balanced gain (i|7|) 
and loss (-i|7|) terms [25|l3T|[32]. 



where 1 < m < N/2 is the position of gain waveguide and rh = N + 1 — m denotes the 
index of the loss waveguide. Although H-pj- is not Hermitian, it has purely real eigenvalues 
and relatively strong VT symmetric phase for a > [33]. We choose random, uniformly 
distributed on-site impurities 7m for 1 < m < N/2; such a choice of random, on-site impu- 
rities leads to a purely real and particle-hole symmetric spectrum [3l|. Motivated by this 



N/2 




(4) 
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result, we compare and contrast the intensity- intensity correlations due to PT-symmetric 
disorder and off-diagonal disorder in tunneling rates. We recall that the position-dependent 
tunneling profile is given by Eq.g, Ca{j) = Co[jiN - j)]'^^^. We introduce the disorder 
in the tunneling rate via Co — Co[l + where the random, position-dependent change 
6{j) is drawn from a uniform distribution with zero mean and variance cr^ and confine our- 
selves to a weak disorder which ensures that the scale-factor 1 + 6{j) is always positive. The 
disorder- averaged, classical, steady-state correlation matrix is defined as |i30i| 



{Iij,t)Iik,t)) 



where I{j,t) is the intensity profile which depends upon the initial state \tpg{0)), and (■ ■ ■ ) 
indicates averaging over different disorder realizations and different relative phases 6 between 
the light input into the two waveguides \p) and \q). Note that since the PT-symmetric 
Hamiltonian, Eq.Q, is not Hermitian, the corresponding time-evolution operator is not 
unitary, and therefore the total intensity Y2j ^Hi ^) ni^-Y not be conserved. For the purpose 
of illustration, we show the results for a uniform (a = 0) = 20 waveguide array with 
initial state |^0(O)) = {\N/2 - 1) + \N/2))/^. The disorder strength is ha / = 0.02 and 
Nr = 10^. The left-hand column in Fig. |4]show the disorder-averaged, steady-state matrix 
Tjk for PT-symmetric, on-site disorder, panel (a), and tunneling disorder, panel (c). 

Traditionally, however, instead of the entire correlation matrix Tjk with independent 
entries, one considers the correlation function g{Ar) = N^'^ XljLi '^jj+Ar = 5'(~^r) with 
N independent entries. This correlation function is able to distinguish between on-site 
potential disorder and off-diagonal tunneling-rate disorder pOj. The right-hand column 
in Fig. |4] shows correlation functions extracted from the steady-state matrix Tj^ for on- 
site PT-symmetric disorder, panel (b), and off-diagonal tunneling rate disorder, panel (d). 
The similarity between the two correlation functions, and their stark contrast with the 
corresponding correlation function for a Hermitian on-site disorder [30], shows that the 
particle-hole symmetry in the spectrum of a disordered system, rather than the origin of the 
disorder, is instrumental in determining the correlation function properties. 
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FIG. 4. (Color Online) Left-hand column shows the disorder-averaged, steady-state, classical 
correlation matrix Fj^ for a uniform array with = 20 waveguides and a weak disorder Ha / — 
0.02; the results are averaged over Nr = 10^ disorder realizations. Panel (a) shows the matrix 
for on-site, PT-symmetric, non-Hermitian disorder; panel (b) shows the matrix for tunneling-rate, 
Hermitian disorder. The right-hand column shows the correlation functions g{Ar) extracted from 
the steady-state, classical correlation matrix for on-site, PT-symmetric disorder, panel (c), and 
off-diagonal, tunneling disorder, panel (d). The similarity between the two results shows that 
the particle-hole symmetry of the disordered energy spectrum is instrumental to the correlation 
function properties. 



V. DISCUSSION 



In this paper, we have explored the effects of disorder in waveguide arrays with non- 
uniform tunneling and non-Hermitian, on-site, PT-symmetric disorder, by focusing on the 
behavior of disorder-averaged, steady-state intensity profile and intensity-intensity correla- 
tions. 

Broadly, we found that the intensity profile I{j, t) is acutely sensitive to the tunneling 
function C^, the initial input state lipeiS^)), its width and its proximity to one of the edges, 
and does not necessarily result in an exponential localization profile that is well known in 
uniform waveguide arrays [9J. In particular, we found that when a > 1 and input state 
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is localized to waveguide mo, a weak disorder suppresses the wave packet reconstruction 
and causes localization of wave packet in its mirror-symmetric waveguide + 1 — mo, as 
well as its initial waveguide mo- We also found that the steady-state intensity encodes 
the phase-information of the initial state lipelO)) for a weak disorder. In all strong 
disorder ha/ Aa^ ^ 1, results in a localized intensity profile that is virtually identical to the 
initial intensity profile. 

Although the localization intensity profile is insensitive to the origin of the disorder, 
on-site or off-diagonal, higher-order intensity correlations depend upon it; in particular, 
Hermitian on-site and off-diagonal disorders lead to qualitatively different intensity correla- 
tion function g{Ar) [30j. Here, we have shown that a Hermitian disorder in the tunneling 
rate, and a non-Hermitian, PT-symmetric, on-site disorder result in nearly identical corre- 
lation functions. Thus, the behavior of the correlation function g{Ar) can be traced to the 
presence or absence of particle-hole symmetry in the spectrum of a disordered Hamiltonian. 

In this paper, we have ignored the effect of Kerr nonlinearity that becomes relevant at 
high electric field amplitudes, and that affects the disorder-induced localization of positive 
and negative energy states differently fSl]. Even in the absence of this nonlinearity, our 
results show that the properties of localization in tunable waveguide arrays, including the 
dependence of saturation time Tg beyond which a steady-state intensity profile emerges on 
the array parameters, are barely explored. Finally, we find that the localization of light in a 
finite lattice can have a complicated dependence on a, the array size and the disorder, and 
further work is necessary to identify the interplay among these factors in determining the 
behavior of light in a coupled array of waveguides. 
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